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Abstract 

We use perturbation theory and the relativistic Cowling approximation to numerically compute 
characteristic oscillation modes of rapidly rotating relativistic stars which consist of a perfect fluid 
obeying a polytropic equation of state. We present a code that allows the computation of modes 
of arbitrary order. We focus here on the overall distribution of frequencies. As expected, we 
find an infinite pressure mode spectrum extending to infinite frequency. In addition we obtain 
an infinite number of inertial mode solutions confined to a finite, well-defined frequency range 
which depends on the compactness and the rotation frequency of the star. For non-axisymmetric 
modes we observe how this range is shifted with respect to the axisymmetric ones, moving towards 
negative frequencies and thus making all m > 2 modes unstable. We discuss whether our results 
indicate that the star’s spectrum must have a continuous part, as opposed to simply containing an 
infinite number of discrete modes. 
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I. INTRODUCTION 


Helioseismology and asteroseismology have provided a wealth of information abont onr 


snn, early type variable stars and even white dwarfs [seeQ, chap. 2], However, for nentron 
stars, this held is still at its beginning, since it is very difhcnlt to measnre their oscillation 
modes in the electromagnetic spectrnm. The gravitational wave spectrnm, thongh, will likely 


allow ns to stndy nentron star oscillations and obtain in 
snch as mas^ radius, rotation rate and equation of state 
Andersson 


brmation about their parameters. 


a. 


In this context, a discovery by 


3j caused some excitement: the so-called r-modes may be unstable at any non¬ 
zero rotation rate of a compact star. Once excited — e.g. during the birth of a neutron star, 
via an accretion process, or through tidal interactions with a bound compact object — it 
may grow to sufficient strength to be detected by an earth-bound gravitational wave detector 
jd]. This mechanism, the so-called CFS instability j^, may cause the fundamental pressure 
mode (/-mode) of a neutron star to become unstable as well. However, this occurs only 
when the rotation frequency of the star is close to the Kepler limit . Many gravity modes 
(gf-modes) can become unstable as well, although their growth is most likely suppressed 
by viscous dissipation a. The .attee are related to either eonrposit.orr gradierrt or hn.te 
temperature, which are not considered in this work. 

In this study we concentrate on the overall distribution of eigenfrequencies rather than 
looking at individual modes. The identihcation of individual modes will be the focus of a 
subsequenct paper. Studying the spectrum as a whole is useful to shed light on questions 
such as: What are the frequency ranges where the star may oscillate at all? How do they 
depend on the physical parameters of the star, such as its rotation rate? Which of them 
may become unstable? Does the spectrum have a continuous part, and are there discrete 
modes within the frequency range of a continuous spectrum? 

We concentrate here on isentropic models, where a star’s equilibrium as well as its per¬ 


turbations can be c 
constant entropy 


escribed with the same one-parameter equation of state p = p{e) and 
BJ (also called barotropic 9|). Deviations from an isentropic model be¬ 
come important only if the angular spin frequency is comparable to or smaller than the 
Brunt-Vaisala frequency, which for neutron stars is of the order of lOOHz |^. 

r-modes have originally been dehned for Newtonian models as inertial modes of axial type. 
For relativistic models this dehnition is not as straightforward any more since the oscillation 
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modes (except dipole modes) consist of a mixture of polar and axial contributions, leading 
to what was called ‘hybrid‘ modes Still, even in the relativistic case with rotation, it 
is possible to assign an overall parity to a mode by looking at the dominant term in an 
inhnite series representation over the harmonic index 1. Lockitch & Friedman js] call these 
modes axial-led or polar-led. If the spherical harmonic with lowest ^ that contributes to the 
velocity field is axial, then the mode was refered to as ’generalized r-mode’ [^. Following 


12l | , this terminology has been made obsolete by a better understanding of the problem. In 


the following we will therefore always use the term “inertial modes”. 

In the relativistic slow rotation approximation 10|, ll3(], their frequencies are found to 
be up to 15% lower than those of their Newtonian counterparts. Using the relativistic 
’traditional’ approximation a similar shift was found [^, but it applies to the frequencies 
in the corotating frame, rather than those in the asymptotic inertial frame. In all cases, the 
shift is smaller for higher order modes. 

In non-rotating stars, inertial modes are degenerate at zero frequency due to the absence 
of their restoring force, the Coriolis force. Rotation breaks this degeneracy, and for isentropic 
stars there is an infinite number of inertial modes confined to a finite frequency range 151. 
The range of frequencies they cover has been shown for Newtonian incompressible stars to 
extend from (—2 — m)i> to (2 — m)z/, where v is the rotational frequency of the star and m 


the azimuthal index 


11| . Later Brink et al. 1^ computed a large set of such modes in the 


same framework and confirmed modes up to 30th order to have frequencies conhned within 

□ 

this range. Ruoff et al. p/7l] studied the inertial mode spectrum for relativistic barotropic 
(as well as non-barotropic) stars in the slow-rotation approximation by including coupling 
of modes up to a maximum harmonic index £max- Next to distinct inertial modes they 
found a continuous spectrum whose width depends on the compactness of the star and on 
£max- They all lie roughly between and uk for m = 0 and between —2h'K and 0 for 
m = 2, where uk ~ 1.915i/ is the breakup frequency of their model. For imax oo the 
authors expect the continuous spectrum to fully cover this range, with individual modes 
still existing inside the continuous spectrum. On the other hand, Lockitch et al. 
their appendix that only individuals modes should be present. Kojima 


9| argue in 


18j suggested that 


the width of the continuous spectrum depends on the range of values that the difference 
between the rotation rate of the star O and the frame dragging uj can take. 

All stars show pressure driven modes (p-modes). The fundamental p-mode, having the 
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lowest frequency, is called the /-mode. They were studied hrst in non-rotating models. 
In Newtonian gravity the perturbation equation describing their frequencies is of Sturm- 
Liouville type [see eg. Il9|, the solution is an inhnite set of modes with frequency range 
unbound from above. General relativity and rotation do not change this general picture, 
and individual frequencies are only slightly affected for a spin frequency up to about half 
the break-up frequency of the star 2^. The /-mode typically has a frequency around 2kHz, 
the lowest order p-mode a somewhat higher one. For a review regarding fluid oscillations of 
neutron stars see 


22[ | and j2]|, [23(] for rotating neutron stars, and [2^ for r-modes. 


Many studies of neutron star oscillations use the so-called Cowling approximation 19| . 


which in its relativistic version allows to work only with fluid displacements and neglect 


metric perturbations 25| . Although p-modes can be computed quite accurately in this 
way 2^, Finn 27| questioned whether neglecting all the metric perturbations is a good 
approximation for modes involving large fluid velocities, such as the p-modes. Lockitch et 
al. [l^ extended this argument to the r-modes, which have similar properties, and disputed 
the results of Ruoff et al. [l^ due to this. Yoshida et al. 281 compared the frequencies of the 
fundamental r-mode for slowly rotating relativistic stars with and without the relativistic 
Cowling approximation and found the difference to be only a few percent of the rotational 
frequency of the star. The same comparison for rapid rotation is not possible due to the 
lack of results with the full spacetime perturbations, but for Newtonian stars the relative 
difference in the frequency is a few percent 28(]. The last two results indicate that even for 
rapidly rotating relativistic stars the Cowling approximation would give a good estimate of 
the r-mode frequencies. 

On the other hand, using the slow-rotation approximation may be more problematic than 
the Cowling approximation. It seems that going beyond the hrst order in the rotation rate of 
the star, results may differ considerably from those obtained in hrst order 29(]. Therefore, in 


this paper we will discard the slow-rotation limit, performing mode calculations for rapidly 
rotating stars using linear perturbation theory in the Cowling approximation in relativistic 
rather than Newtonian gravity. Given the results discussed above, we expect to hnd an 
inhnite number of pressure modes with frequencies on the order of kHz and higher, and 
possibly an inhnite number of inertial modes with frequencies dependent on the rotation 
rate of the star. 

In section [TTl we formulate the problem and discuss the numerical setup. In section IHII 
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we investigate axisymmetric perturbations, first in the non-rotating case for comparison 
with known results, and then for the rotating case. Results for non-axisymmetric perturba¬ 
tions and a comparison with the axisymmetric ones follow in section IIVI concluding with a 
discussion of the results in Section IVl 


II. PROBLEM SET-UP 


A. Equilibrium background 


Neutron stars have been observed to spin up to millisecond periods. This results in 
signihcant deviations from spherical symmetry both for the fluid conhguration of the star 
and for the spacetime. Equilibrium solutions will thus not be spherically symmetric, but 
axisymmetric at most. The general form of an axisymmetric metric describing a rotating 
body can be written as: 


ds‘^ = + e^“ {dr'^ + r‘^d6‘^) + sin^ 9 {dcj) — cudtY ■ 


( 1 ) 


where r,6,(j) are quasi-isotropic coordinates, reducing to isotropic ones if the rotation rate 
goes to zero[see l23(]. Restricting ourselves to uniform rotation with frequency U, the 4- 
velocity of the fluid is given by u" = 0, 0, U}, with an energy-moment tensor = 

P9fiu + (p + e) UfJJy, assuming that the star consists of an ideal fluid. For the equation of 
state we use a simple polytropic model of the form p = k x p^. Using a realistic equation 
of state in tabulated form would not affect our analysis or our numerical procedure. In this 
paper, we restrict ourselves to a polytropic equation of state to facilitate comparison of our 
results with previous studies by other authors. 

We use the RNS code of Stergioulas & Friedma n |30l | for computing the background 
models. The same code was also used by Font et al. j^]; we will occasionally refer to their 
results for comparison. Table [I] summarizes some of their models which we will be using 
here. For all models, the parameters in the equation of state are k = 217.856km^ and 7 = 2; 
the central density is always Cc = 0.894 x lO^^g/cm^. 
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TABLE I: The models of Font et al. 


20( 1 which we refer to in this paper. Next to the resulting 


mass and radius we show both the angular frequency, that is used in our analysis, and the natural 
rotation rate. 


Model 

Gravitational Mass 

Radius 

n 

Rot. rate (u) 

BUG 

1.4Mq 

14.15 km 

0 

0 

BUI 

1.432M0 

14.51 km 

2.185 kHz 

348Hz 

BUG 

1.627M0 

17.25 km 

4.984 kHz 

793Hz 


B. Perturbations 


Assuming small deviations for the fluid variables we study linearized perturbations on 
these stationary conhgurations. Since the background is not spherically symmetric it is not 
helpful to decompose the perturbations into spherical harmonics. Instead we exploit the 
axisymmetry of the background by writing the perturbation quantities as 

(2a) 

1 


5Un = 


p + e 


/i, /2, /sje 


im(f> 


(2b) 


where H,fi are functions of t, r and 6. The reason for this particular choise will become 
clear in the next subection fill CD . The perturbation of the energy-density is related to the 
pressure perturbation through 


Sp=^Se = C^Se, 
de 


where Cg is the speed of sound. 


As mentioned earlier we use the Cowling approximation, ( 2^ after 19|, iM]) neglecting 
the metric perturbations. This will not allow us to calculate any damping of the modes due 
to emission of gravitational waves, but we can estimate the oscillation frequencies and study 
the overall structure of the spectrum (see also section [I]). 

The equations that describe the behavior of the perturbed quantities arise from the 
perturbed form of the conservation of energy-moment (equations of motion for the fluid): 


h - T^.hT,,) = 0 


(3) 
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In general, these yield four independent partial differential equations which are of first order 
in time and space. Time-evolving sets of lengthy equations can often lead to instabilities 


[see eg. 1^. We perform mode calculations by assuming a harmonic time-dependence 
for all four variables (e.g. H{t,r,6) = H{r,6)e^‘^^), searching for frequencies which allow 
non-trivial solutions of the perturbation equations. 


C. Numerical procedure 


Here we briefly describe the numerical method used in this work. We discretize the system 
of equations at every grid point, including the boundaries by making use of the boundary 
conditions. This results in a system of linear equations of the form 

A ■ X = 0 (4) 


where the unknowns X are the discrete perturbation quantities at all grid points, and the 
coefficient matrix A represents the equations resulting from the perturbation equations at all 
grid points and from the boundary conditions. A is a highly sparse matrix. The components 
of A depend on the frequency parameter a. It is now possible to examine the condition of 
A, searching for values of a which make A degenerate, since this is the only way that Eq. 
()T|) allows non-trivial solutions for the perturbation quantities. 

In the case we study here, we may rewrite Eq. (jl]) as an eigenvalue problem of the form 

A-X = iaX (5) 


Note that this is a special case and not generally possible; it will probably not work for the 
general perturbation equations which result if one does not use the Cowling approximation. 

We may now use a standard routine for finding eigenvalues of Eq. ([5]). We use the 
routine cg.f from the EISPACK package of the NETLIB libraries which handles general 
complex matrices by use of the QR-algorithm. The characteristic frequencies of the star’s 
oscillation modes are obtained as eigenvalues of Eq. ([5]), the perturbation quantities X are 
the corresponding eigenfunctions. It has turned out that a two-sided differencing scheme 
does not give consistent results (see |33|). We therefore use a one-sided differencing scheme 
for all hrst derivatives throughout this work. 

We need to take into account boundary conditions at the center and at the surface of 
the star, and a regularity condition in the angular direction. At the center, all variables are 
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required to vanish by the regularity condition there. We implement this condition simply by 
setting all variables to zero at r = 0. At the surface, all perturbation variables must vanish: 
this follows from the dehnition of the variables given in Eqs. ([2]), since p = e = drP = 0 (for 
polytropic equations of state) and Sp ~ drP (to hrst order in a Taylor expansion). For the 
special case of m = 0, the pressure perturbation is not required to vanish at the center; its 
value must then be determined directly through Eqs. ([3]). 

For the boundaries in the 6'-direction one may use the fact that the rotational axis itself 


is special: 
conditions 


or m > 0, all variables have to vanish on the rotational axis due to regularity 


371]. This can be used directly to set values for the discretized variables there. 


For m = 0 we may construct a grid in such a way that no grid points fall on the rotational 
axis. In this case, one must use the symmmetry condition to construct a boundary condition 
in the 6'-direction. This technique is also applicable for m > 0, keeping in mind that the 
symmetry is different for even and odd values of m. We have used both techniques for the 
m = 2 case and found no signihcant difference. In order to make the code applicable both 
in the axisymmetric and the non-axisymmetric case, we implemented the second possibility. 

The RNS code itself has hnite resolution and thus introduces some numerical error into 
the results. We will always see a combination of numerical error coming from the RNS code 
and from our own code. When we attempt to study convergence properties we therefore use 
a high but hxed resolution of the RNS code, varying only the resolution of our eigenvalue 
code. In order to avoid the use of interpolation and the additional error associated with it, 
we run our code only with grid spacings which are multiples of the hxed grid spacing we 
have chosen for the RNS code. 

III. AXISYMMETRIC PERTURBATIONS 


We hrst study axisymmetric perturbations as a test for our code. We use the correspond¬ 
ing equations in their general form, including all terms with arbitrary m. Setting m = 0 we 
select perturbations that are symmetric around the rotation axis, including all £ > 0 contri¬ 
butions. Axisymmetric perturbations have been studied by Font et al. 20j using nonlinear 
time-evolution; we will turn to their results for comparison. In order to gain a better under¬ 
standing of the underlying physics and numerics, we hrst solve Eqs. ([3]) for the non-rotating 
case where expressions may be reduced considerably. We then proceed to include rotation 
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at arbitrary rotation rate. 


A. No rotation 


Setting = 0 in Eqs. ([3D and using harmonic time-dependence the system of equations 
for the variables dehned in Eqs. [2] takes the form: 


, ^ c; f 1 (dh mse \ a/i 

r^sin^^ ^ \ r‘^ \ d9 sin^ dr 


Sd'j Idp 2 

- - H- - + - 

2 dr 2 dr r 


fi 


tafs = lafi = 


1 dH 


iaf2 = 


1 dH 


( 6 ) 


UOdr' d9 

In this simple case, one variable {H) would be sufficient to describe the whole oscillation 
problem since one can transform the system of Eqs. (ED into a single second order equation. 
Solving the latter gave the same results as solving Eqs. (ED for = 0, so we will always 
present results from the full hrst-order system for consistency. 

Table mi shows the lowest resulting eigenfrequencies for the model BUO with increasing 
resolution in the radial direction, starting from 21 points. Increasing the resolution in the 
9 direction does not seem to affect the results, and it is computationally expensive. We 
therefore keep the ^-resolution hxed for the results shown in Table [TTl From now on we 
will always use the physical oscillation mode frequency / rather than the angular frequency 
(a = 271 f), unless stated otherwise. 
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FIG. 1: The computed eigenfrequencies (here the i = 2 /-mode for BUO) follow an inverse power 
law. We may extrapolate to ^ oo to obtain the frequencies for nominally infinite radial 
resolution. 

Figure [1] shows a plot of the (/ = 2, m = 0) /-mode frequencies from Table HI] as a function 
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TABLE II: The frequencies (^) in Hz of the first three fundamental modes (/) and the lowest 
pressure mode (p^) for axisymmetric oscillations of BUG, for several resolutions in the radial direc¬ 
tion and a fixed number of points in the angular direction. In the next to last row the extrapolated 
values for ’infinite’ radial resolution (see figured]) are listed, while the last row shows the results 


of Font et al. 


201 ]. 


ne Tir 

fe=i 

fl=2 

fe=o 

P\=2 

10 21 

1377 

1974 

3076 

4627 

26 

1367 

1943 

3032 

4551 

51 

1344 

1872 

2917 

4344 

101 

1328 

1830 

2840 

4197 

oo 

1317±2 

1794±3 

2787±11 

4102±26 

Font et al. [20] 

1335 

1846 

2706 

4100 


of the radial resolution. They follow an inverse power law of the form / = /oo + ^//'^r + 
..., exhibiting hrst order convergence as one expects for one-sided differences. We may 
thus extrapolate the computed values to obtain an eigenfrequency foo at nominally inhnite 
resolution. In general, the relative mean deviation of the htted numbers is less than 1%. 
Note that there is an additional systematic error of the same general magnitude, resulting 
from the hnite accuracy of the background quantities. 

Since the background is now spherically symmetric, we can assign a dehnite value of £ to 
the oscillation modes. Comparing the extrapolated frequencies in Table [TTl with results by 
Font et al. , we conclude that they correspond to the fundamental modes that correspond 
to £ = 0,1, 2 and the hrst pressure mode for i = 2, with an agreement up to a few percent. 
The remaining differences are likely due to the hnite resolution of the background calculation. 

The method remains basically the same when we turn on rotation. We thus expect to 
obtain the huid modes with similar accuracy even for rapid rotation. The picture may 
change, though, for the rot at ion-driven modes: These are degenerate at zero frequency in 
the non-rotating case. Therefore, there are no results for the non-rotating case that can be 
used as an indication for the accuracy of our computational results. 
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B. Rapid rotation 


With rotation being included the equations of motion for the fluid become quite lengthy. 
Schematically they take the form 


n 


dH 

idfs 


— \ H, /s, /i, /2, 


dfi df2 
dr ’ 89 


,-AC: 


dt 




dt 


- C (f f 
^ = 1 : (f f 


where A,B,Z,n are functions of the background quantities (provided in Appendix!^. How¬ 
ever, the frequency a now appears in the matrix C on the right hand side. Therefore, 
this form is not suitable for performing an explicit eigenvalue calculation as described in 
Subsection III Cl We thus define new variables: 


F = UH + AClh 
V = zclh + bh, 


which together with /i, /2 form the set of variables we will use from now on. The full set 
of equations is provided in Appendix 

In Fig. [ 2 ] we show all solutions returned by the eigenvalue code for axisymmetric pertur¬ 
bations of model BUI, using a low radial resolution of = 25. As far as results for pressure 


TABLE III: Frequencies (in kHz) of two fundamental m = 0 
polytropic models BUI and BU6, compared with the results of 


pressure-driven oscillations for the 


20 ( 1 . In the non-rotating limit these 


correspond to I = 0 and I = 2 modes. Again, our values have been extrapolated to nominally 
infinite radial resolution. Convergence in this case was not as clean as shown in Fig. (H especially 
for the rapidly rotating model BUG. This is reflected in larger uncertainty estimates. 



fi=o (BUI) 

/i=2(BUl) 

//=o(BU6) 

/z=2(BU6) 

This paper 

2.720±20 

1.834±25 

2.292±193 

1.718±85 

Font et al. [20] 

2.657 

1.855 

2.456 

1.762 
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Inertial 



/[kHz] 


FIG. 2: A histogram of the number (per frequency bin) of all solutions returned by the eigenvalue 
code for the m = 0 eigenvalue problem of the model BUI. The frequency bins have a constant 
width of 10 Hz for the inertial mode frequencies and 1 kHz for the pressure mode frequencies. Note 
the change in frequency scale at IkHz. 


modes are available from the literature we list them along with our results in Table IIIII No 
previous results are available for axisymmetric inertial modes of rapidly rotating relativistic 
stars; see Ruoff et al. ITj] for modes of slowly rotating stars. Very recently Dimmelmeier et 


al. 


34( 1 published frequencies of the three strongest axisymmetric inertial modes in the con¬ 


formally flat approximation. They all lie inside the corresponding inertial mode spectrum 
(see below). 

We notice from Fig. [2] that the eigenvalues are clustered into two groups, one above 
lOOOHz and one below. These correspond to the expected frequency ranges for pressure 
modes and inertial modes. The latter range is more densely populated; for example, we see 
more solutions between, say, 500 and 520 Hz than between 2 and 5 kHz. This is not as 
surprising as it may seem at hrst: according to theory (see eg. for an overview table), 
there is an infinite number of pressure modes, with frequencies extending to infinity. For 
the inertial modes, though, one expects an inhnite number of modes as well, but conhned to 
a well defined frequency range. According to our computation this range appears to extend 
from 0 to about 600 Hz for the model BUI. 

In Figs. [3] and m we show the p—mode and inertial mode ranges separately for increas¬ 
ing radial resolution. In both ranges the number of frequency eigenvalues increases with 
increasing resolution. There is an important distinction, however: in the p-mode range. 
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FIG. 3: A histogram of the number of p-modes 
per frequency bin computed for the m = 0 
eigenvalue problem of the model BU1. The bins 
are equally sized at lOkHz 



FIG. 4: Same as Fig. O but for the inertial 
modes of the same problem. The bin size is 
now 25Hz 


higher frequency ranges are increasingly populated as the resolution increases, while the 
population of low frequencies approaches a limiting value. We observe that as the radial 
resolution rij. increases, the number of computed eigenvalues increases as 10 x u,,. The lowest 
p-mode frequency, which belongs to the m = 2 /-mode, is roughly the same for all resolu¬ 
tions. It shifts in frequency by a few hundred Hz, which is just a few percent of the bin 
size. Typically the number of solutions per frequency bin approaches a limit value a like 
Win ^ a + P/ur with some constant f3. 

For each stellar model the maximum value of frequency eigenvalues tends to infinity for 
increasing resolution. This is similar to the situation for quasi-normal modes of black holes, 
with an inhnite number of modes which is not conhned to a hnite part of the complex 
frequency plane (see 0 ). 

In the inertial mode range, on the other hand, the frequency range does not change as the 
resolution increases; in fact, the upper limit (600Hz for model BUI) is quite robust. Instead, 
the population increases fairly homogeneously over the whole frequency range. The total 
number of solutions calculated develops, just as for the p-modes, like 10 x n^. The number 
of points per bin grows linearly; for the bin e.g. around 400Hz, as Win — 2 x n^. This is 
a strong indication that an infinite number of solutions exists in this frequency range. It 
would seem likely that there is an infinite number of physical modes in this range as well. 
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While the upper limit for the frequency range does not depend on the resolution used 
in the numerical calculation, it actually depends linearly on the rotational frequency of the 
star, as shown in Fig. O The linear fit reveals Cmax = 1-674 x hi (where a and hi are now 
the angular frequencies), with a negligible statistical error. 



FIG. 5: The highest inertial mode angular frequency for m = 0 as a function of the star’s rotational 
angular frequency, starting at the model BUO (no rotation) and moving up to model BUI (2185Hz). 

The picture in Fig.[3]is actually quite common for numerical studies of oscillation spectra: 
As one increases the resolution, frequency ranges at increasingly higher frequency become 
populated, while there is a hnite limit for lower frequency ranges. Usually one expects that 
at low resolution, only modes with low frequencies can be computed reliably, since they are 
the only ones which can be resolved sufficiently well. Numerical solutions at high frequencies 
are likely spurious and cannot be trusted. With increasing resolution, the number of reliable 
solutions increases and their range extends to higher frequencies. 

On the other hand, if we are confronted with a hnite frequency range which shows an 
increasing number of solutions with increasing resolution, such a distinction cannot be made 
in a meaningful way. It is therefore not clear which of the numerical solutions in the frequency 
range corresponding to inertial modes should be considered physical solutions, and which 
should be discarded as numerical artefacts. 

Just as in the non-rotating case, one can establish convergence for any hnite value of 
the frequency in the upper range (Fig. [3l), since there is a hnite number of distinct modes 
in any hnite frequency interval. The limit values very likely correspond to physical modes 
of the star. In the lower frequency range (Fig. 0]) it is not even clear how to establish a 
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correspondence between eigenvalues at increasing resolutions, let alone define convergence 
and establish a correspondence to physical modes. 

However, the number of eigenvalues increases linearly with the number of grid points 
used in the calculation. We therefore find it unlikely that this part of the spectrum contains 
a finite number of discrete oscillation frequencies. If an infinite numbers of eigenvalues is 
confined to a finite frequency range, then they must have at least one accumulation point. 
It is possible that there is only a finite number of accumulation points, even only one. An 
example for such a system is an atom which has an infinite number of bound states confined 
to a finite energy interval, with zero energy as the only accumulation point. However, as 
we increase the resolution of our calculations, the number of eigenvalues increases rather 
evenly throughout the whole frequency range in question, and there is no sign that there 
are parts where the increase might saturate at a finite level. We thus consider it likely that 
the oscillation frequencies form a dense set throughout the inertial mode frequency range. 
The spectrum may also be continuous over this whole range or in parts of it, but this can 
neither be confirmed nor excluded based on the results of our code. 

IV. NON-AXISYMMETRIC PERTURBATIONS 

For m > 0 the picture is more complicated but also more interesting than the axisym- 
metric case since the m = 2 modes are the ones most unstable to gravitational radiation. 
Solving this system for m = 2 and the BUI model in the way described in the previous 
section we find a set of eigenvalues containing both positive and negative frequencies. In the 
axisymmetric case negative eigenvalues are equivalent to the positive frequency solutions. 
Breaking axial symmetry changes frequencies (for an asymptotic observer) towards negative 
values, so that much or even all of the frequency range corresponding to inertial modes 
becomes negative. 

In Newtonian gravity (see e.g. Unno ef a/. [l|), a polar mode of order n, harmonic index 
i, and frequency Uo, splits under rotation to 2^ + 1 modes with frequencies changing to 
(T = (To — mflEne (+0(U^)) where Ene is a function depending on the eigenfunction of each 
mode. For low order pressure modes the value of this function is about 0.1, so one would 
need rotational frequencies close to the Kepler limit for the frequency to change sign. 

However, there are modes which have a < 0 for any rotation rate, such as the r-modes. 
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which, in the Newtonian, slow-rotation limit, have a = This is the pictnre 

we hnd for m = 2, with all inertial modes having negative freqnencies. Snch a change of 
sig n is often used to mark the onset of the instability for the corresponding mode [see eg. 


231 . 


In Fig. [5] we saw that the upper cutoff frequency for the inertial mode range grows 
linearly with the rotational frequency of the star. In the non-axisymmetric case, there 
is also a dependence on m. For the series of models listed in Table [I] we obtain a least 
square fit of cTmax = 12(1.6 — 1.05m), with statistical errors of about 0.1 for both values. 
For less relativistic models with a central energy density 1/10 that of BUI, this changes to 


n 




o'ma.x = 12(1.93 — 1.03m), close to what Lindblom & Ipser [Utj predicted and Brink et al. 
found for Newtonian stars where cXmax ~ f2(2 — m). 

In the axisymmetric case a zero-frequency mode in the rotating frame has 0 Hz for an 
inertial observer as well, corresponding to the mid-point of the inertial mode frequency 
range. For non-axisymmetric perturbations the center of the spectrum is expected to shift 
to —mf2, which indeed appears to be the case in our results. The corresponding inertial mode 
spectrum for a specihc azimuthal index m > 0 is however not completely symmetric since 
modes with different order n or harmonic indices ^ have different frequency changes. This 
change is rather small for pressure modes. As the azimuthal index m can take arbitrarily 
large values, the inertial modes can, according to the above, reach arbitrarily large negative 
frequencies. Since these are equivalent to positive frequencies with a phase difference, the 
frequency range of inertial modes will overlap with that of pressure modes. An overview of 
the spectra for three different values of the azimuthal index m can be seen in Fig. [6l 

The oscillation frequencies are influenced not only by changes in the equilibrium conhg- 
uration of the star, but by other effects as well, such as the frame dragging at the surface 
and the center 36(|. These also scale linearly with 12. Knowing the range of inertial mode 
frequencies as a function of the star’s rotational period in advance may be quite helpful for 
actual observations of gravitational radiation emitted by these oscillations. The frequencies 
of individual modes, such as the fundamental r-mode, are not examined here, they will form 
the focus of a forthcoming paper. 
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FIG. 6: The histogram of both positive and negative solutions of the eigenvalue code for m = 0, 2, 
and 4. The values for the inertial modes (middle) are scaled differently for better presentation. 

V. CONCLUSIONS 

We have presented a code to numerically compute oscillation frequencies of rapidly ro¬ 
tating relativistic stars for both axisymmetric and non-axisymmetric perturbations, using 
the relativistic Cowling approximation. 

For the polytropic equations of state that we have employed here, we found an inhnite set 
of pressure modes with a range of frequencies bounded only from below at about 2kHz. In 
addition, there is a presumably inhnite set of solutions at lower frequencies in a well dehned 
range which corresponds to inertial modes. In the axisymmetric case (m = 0) this frequency 
range is symmetric around zero. It extends to a maximum frequency which depends linearly 
on the star’s rotation rate Vstar and approaches 2ustar for less relativistic stars. For non- 
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axisymmetric perturbations oscillation frequencies change towards negative numbers. This 
affects pressure modes only slightly, while the inertial modes now all have frequencies below 
zero. 

It is not clear whether this low-frequency part of the spectrum is discrete or continuous, 
or a combination of both. The dense and somewhat uneven distribution of the eigenvalues is 
an indication that the actual structure of the physical spectrum may be quite complicated. 
Further investigation is needed into the question how a spectrum with a continuous part can 
be studied numerically. Furthermore, explicit time evolution of perturbations of the same 
conhgurations could provide an independent means of studying their spectrum and checking 
the results we have presented here. Identihcation of individual modes in this part of the 
spectrum will be studied in a forthcoming paper. 

Applying realistic equations of state requires only slight modihcation of the code and will 
shift the frequencies to some extent, but the overall picture should remain the same. The 
effect of a non-barotropic equation of state on the inertial mode spectrum presents a very 
interesting question for further study. 

Finally, for an accurate calculation of the oscillations frequencies, the disposal of the 
Cowling approximation will be necessary. This is important not just to avoid deviation of 
the calculated frequencies introduced by the Cowling approximation, which we expect to be 
small, but also to ensure the consistency of the problem (see Section [T|). The major problem 
lies not in the perturbation equations becoming more complex, but in dehning the boundary 
conditions in the spacetime outside the star. Observation of potentially unstable modes with 
the now-operational gravitational-wave detectors is possible only if precise understanding 
of these modes and their frequency range is available as a basis for the analysis of the 
gravitational wave data. 
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APPENDIX A: THE FULL SYSTEM OF EUATIONS FOR THE FLUID PER¬ 


TURBATIONS 
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dtf2 = 
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and dr, dg stand for the partial derivatives with respect to the radial and angular coordinate. 


II. DISCUSSION OF SOME NUMERICAL DIFFICULTIES 

The Cowling approximation does not include the damping of oscillations by emission of 
gravitational waves, and no other damping mechanism is included in our models. Therefore, 
all eigenfrequencies should be purely real. However, the code we use is designed for calcu¬ 
lating complex frequencies. The imaginary parts of the computed values should therefore 
turn out to be close to zero, on the order of machine accuracy (~< 10“^^). This is indeed 
the case for the majority of the solutions we found. Some of the eigenvalues in the frequency 
range of modes, however, have a non-negligible imaginary part of ~< 1% of the real part. 

We are thus facing the question if some error occurred in the implementation of the per¬ 
turbation equations into the code. The perturbation equations where hrst computed using 
Maple. The resulting expressions were then simplihed by hand and implemented numerically. 
As an alternative, we inserted the Maple output directly into the code without any simplih- 
cation. The spurious imaginary parts of inertial mode frequencies then became larger, rather 
than smaller. We are therefore conhdent that they are not the result of some actual error in 
implementing the equation, but a consequence of truncation error, especially in conjunction 
with a bad numerical condition of the coefficient matrix A. 
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Deviations from the expected convergence were found, while eigenfunctions were also seen 
to be changing form in the rotating case. Testing our code by relabeling the eigenvectors 
revealed the same results, which strengthens its correctness. Still, considerably more work 
is required to clarify the question how numerical results can help to reveal the complicated 
structure of the spectrum we are studying here. 
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